*DECK MOMFLZ
      SUBROUTINE MOMFLZ
C==============================================================MOMFLZ
C     THIS ROUTINE UPDATES RORU, RORV, AND RORW
C     DUE TO CONVECTIVE FLUXING OF MOMENTUM IN THE Z DIRECTION
C----------
C     CALLED BY LAVA
C----------
C     INPUT: ROWU, ROWV, ROWW, DX, DY, DZ, RC, RT
C     OUTPUT: RORU, RORV, RORW
C==============================================================MOMFLZ
      INCLUDE 'COML.h'
C----------
      DOUBLE PRECISION
     1     RDZCT,WC,SS,COUR,WTC,WNTR,WNTF,DELT,RONC,RONB,RONT,
     2     ROUC,ROUT,ROVC,ROVT,ROWB,ROWT,RTF,RTR,SSB,SST,DRB,DRT
C
C==============================================================MOMFLZ
      IF(NX.EQ.1) GO TO 300
      DO 210 IC=IC1W-1,IC2
      WNTR=(DX(IC)*WN(IC+1)+DX(IC+1)*WN(IC))*HRDXCR(IC)
      SS=SIGN(ONE,WNTR)
      RDZCT=TWO*HRDZCT(IC)
      DELT=DZ(IC+NXYT)*RDZCT
      COUR=WNTR*DT(IC)*RDZCT
      WTC=HALF*(DELT+ADC*(ONE-DELT+SS)+BDC*COUR)
      RONT=(DX(IC)*RON(IC+NXYT+1)+DX(IC+1)*RON(IC+NXYT))*HRDXCR(IC)
      RONC=(DX(IC)*RON(IC+1)+DX(IC+1)*RON(IC))*HRDXCR(IC)
      ROUT=RONT*UN(IC+NXYT)
      ROUC=RONC*UN(IC)
      RTR=(DX(IC)*RT(IC+1)+DX(IC+1)*RT(IC))*HRDXCR(IC)
      RORWU(IC)=(ROUC*WTC+ROUT*(ONE-WTC))*RTR*WNTR
  210 CONTINUE
      DO 240 IC=IC1U,IC2
      RORU(IC)=RORU(IC)-DT(IC)*(RORWU(IC)-RORWU(IC-NXYT))*RDZ(IC)
  240 CONTINUE
  300 CONTINUE
C==============================================================MOMFLZ
      IF(NY.EQ.1) GO TO 200
      DO 110 IC=IC1W-NXT,IC2
      WNTF=(DY(IC)*WN(IC+NXT)+DY(IC+NXT)*WN(IC))*HRDYCF(IC)
      SS=SIGN(ONE,WNTF)
      RDZCT=TWO*HRDZCT(IC)
      DELT=DZ(IC+NXYT)*RDZCT
      COUR=WNTF*DT(IC)*RDZCT
      WTC=HALF*(DELT+ADC*(ONE-DELT+SS)+BDC*COUR)
      RONT=(DY(IC)*RON(IC+NXYT+NXT)+DY(IC+NXT)*RON(IC+NXYT))*HRDYCF(IC)
      RONC=(DY(IC)*RON(IC+NXT)+DY(IC+NXT)*RON(IC))*HRDYCF(IC)
      ROVT=RONT*VN(IC+NXYT)
      ROVC=RONC*VN(IC)
      RTF=(DY(IC)*RT(IC+NXT)+DY(IC+NXT)*RT(IC))*HRDYCF(IC)
      RORWV(IC)=(ROVC*WTC+ROVT*(ONE-WTC))*RTF*WNTF
  110 CONTINUE
      DO 140 IC=IC1V,IC2
      RORV(IC)=RORV(IC)-DT(IC)*(RORWV(IC)-RORWV(IC-NXYT))*RDZ(IC)
  140 CONTINUE
  200 CONTINUE
C==============================================================MOMFLZ
      DO 10 IC=IC1,IC2
      WC=HALF*(WN(IC)+WN(IC-NXYT))
      SS=SIGN(ONE,WC)
      COUR=WC*DT(IC)*RDZ(IC)
      WTC=HALF*(ONE+ADC*SS+BDC*COUR)
      RONB=(DZ(IC)*RON(IC-NXYT)+DZ(IC-NXYT)*RON(IC))*HRDZCT(IC-NXYT)
      RONT=(DZ(IC)*RON(IC+NXYT)+DZ(IC+NXYT)*RON(IC))*HRDZCT(IC)
      ROWB=RONB*WN(IC-NXYT)
      ROWT=RONT*WN(IC)
      RORWW(IC)=(ROWB*WTC+ROWT*(ONE-WTC))*RC(IC)*WC
   10 CONTINUE
C---------- EXTRAPOLATE OUTFLOW BOUNDARY ----------------------MOMFLZ
      DO 20 ICBT=1,NBT
      ICB=ICBT
      ICT=ICBT+NXYZT-NXYT
      SSB=SIGN(ONE,W(ICB))
      SST=SIGN(ONE,-W(ICT-NXYT))
      DRB=(DZ(ICB)+DZ(ICB+NXYT))/(DZ(ICB+NXYT)+DZ(ICB+NXYT+NXYT))
      DRT=(DZ(ICT)+DZ(ICT-NXYT))/(DZ(ICT-NXYT)+DZ(ICT-NXYT-NXYT))
      RORWW(ICB)=RORWW(ICB+NXYT)+HALF*(ONE-SSB)*DRB
     1           *(RORWW(ICB+NXYT)-RORWW(ICB+NXYT+NXYT))
      RORWW(ICT)=RORWW(ICT-NXYT)+HALF*(ONE-SST)*DRT
     1           *(RORWW(ICT-NXYT)-RORWW(ICT-NXYT-NXYT))
   20 CONTINUE
C----------
      DO 40 IC=IC1W,IC2
      RORW(IC)=RORW(IC)-DT(IC)*(RORWW(IC+NXYT)-RORWW(IC))*TWO*HRDZCT(IC)
   40 CONTINUE
C==============================================================MOMFLZ
      RETURN
      END
